package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T62;
import scu.maqiang.fes.Materials;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Elasticity2DT62 {

	public static void main(String[] args) {
		Mesh2T6 mesh = new Mesh2T6().square(80, 10, T3Type.Left);
		mesh.scale(8.0, 1.0);
		double E = 1.0e10;
		double Nu = 0.3;
		double[][] D = Materials.ElasticConstituteMatrix(E, Nu, BVPType.PLANESTRESS);
		FES2T62 fs = new FES2T62(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[][][] {D}, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		fs.assembleSource(new double[] {0.0, -1.0e3}, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		fs.applyBC_MBN(A, Direct.All, 4);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		double[] x = new double[fs.GetNdof()];
		solver.PCGSSOR(RHS, x, 1.5, 1);
		double[][] uv = new double[2][mesh.getNv()];
		fs.extractComponent(x, uv);
		System.out.println("Solution: ");
		System.out.println("u min: " + MVO.min(uv[0]) + "\t" + "max: " + MVO.max(uv[0]));
		System.out.println("v min: " + MVO.min(uv[1]) + "\t" + "max: " + MVO.max(uv[1]));
		mesh.toTecplot("Displacement.dat", uv);
		
		double[][] FD = new double[4][mesh.getNv()];
		double[][] SD = new double[6][mesh.getNv()];
		fs.computeFirstAndSecondDerivative(uv, FD, SD);
		System.out.println("First Derivative: ");
		System.out.println("ux  min: " + MVO.min(FD[0]) + "\tmax: " + MVO.max(FD[0]));
		System.out.println("uy  min: " + MVO.min(FD[1]) + "\tmax: " + MVO.max(FD[1]));
		System.out.println("vx  min: " + MVO.min(FD[2]) + "\tmax: " + MVO.max(FD[2]));
		System.out.println("vy  min: " + MVO.min(FD[3]) + "\tmax: " + MVO.max(FD[3]));
		System.out.println("Second Derivative: ");
		System.out.println("uxx  min: " + MVO.min(SD[0]) + "\tmax: " + MVO.max(SD[0]));
		System.out.println("uyy  min: " + MVO.min(SD[1]) + "\tmax: " + MVO.max(SD[1]));
		System.out.println("uxy  min: " + MVO.min(SD[2]) + "\tmax: " + MVO.max(SD[2]));
		System.out.println("vxx  min: " + MVO.min(SD[3]) + "\tmax: " + MVO.max(SD[3]));
		System.out.println("vyy  min: " + MVO.min(SD[4]) + "\tmax: " + MVO.max(SD[4]));
		System.out.println("vxy  min: " + MVO.min(SD[5]) + "\tmax: " + MVO.max(SD[5]));

		Mesh2T3 mesht3 = new Mesh2T3().square2D(640, 80, T3Type.Left);
		mesht3.scale(8.0, 1.0);
		Elasticity2DT32 e2dt32 = new Elasticity2DT32(mesht3);
		e2dt32.elastoStaticProblem(BVPType.PLANESTRESS, new double[] {1.0e10, 0.3}, new double[] {0.0, -1.0e3}, 
				new int[][] {new int[] {4}}, new Direct[] {Direct.All}, new double[][] {new double[] {0.0, 0.0}}, 
				null, null, 
				null, null, null, "MyPlaneStress2DT323.dat");
//		Mesh2T6Generator mg = new Mesh2T6Generator();
//		Mesh2T6 mesh = mg.square(80, 10, T3Type.Left);
//		mesh.scale(8.0, 1.0);
//		double E = 1.0e10;
//		double Nu = 0.3;
//		double[][] D = Materials.ElasticConstituteMatrix(E, Nu, BVPType.PLANESTRESS);
//		FES2T62 fs = new FES2T62(mesh);
//		SRMatrix A = new SRMatrix(fs.GetNdof());
//		double[] RHS = new double[fs.GetNdof()];
//		fs.assembleElasticStiff(new double[][][] {D}, BVPType.CONSTITUTE_MATRIX_COMMON, A);
//		fs.assembleVolumeForce(new double[] {0.0, -1.0e3}, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
//		fs.applyBC_MBN(A, Direct.All, 4);
//		fs.applyBC_MBN(RHS, Direct.All, 0.0, 4);
//		NewIterSSolver solver = new NewIterSSolver(A);
//		double[] x = new double[fs.GetNdof()];
//		solver.PCGSSOR(RHS, x, 1.5, 2);
//		double[][] uv = new double[2][mesh.getNv()];
//		fs.extractComponent(x, uv);
//		System.out.println("Solution: ");
//		System.out.println("u min: " + MVO.min(uv[0]) + "\t" + "max: " + MVO.max(uv[0]));
//		System.out.println("v min: " + MVO.min(uv[1]) + "\t" + "max: " + MVO.max(uv[1]));
//		mesh.toTecplot("Displacement.dat", uv);
	}

}
